	PROGRAM SCRUGGS
! program to compute equilibria for Pallage-Scruggs-Zimmermann UI measurement paper
! this version: 20.5.2007 (CZ)
	REAL*8 MAXAST,INCAST,BETA,CONVERGE,PI,MPI,VTEST,BAD,FTEST,INT,MR
	INTEGER I,J,K,L,M,N,O,P,Q,R,NASTS1,NASTS2,ITER,ITERV,ITERF,J1,K1,L1,M1,N1,O1,YY
	REAL*8 HBAR,SIGMA,MHBAR,MSIGMA,MBAD,TUTU
	REAL*8 THETA1,TAU1,MTAU1,THETA2,TAU2,MTAU2,IS,TOTO,RHO,MRHO,TOTO1,TOTO2
	PARAMETER (BETA=0.999165D0,MAXAST=4.0D0,INCAST=0.025D0,CONVERGE=0.001D0)
	PARAMETER (YY=32,INT=0.0D0,MR=1.0D0+INT)
	PARAMETER (BAD=-100000.0D0,HBAR=0.45D0,SIGMA=0.67D0)
	PARAMETER (NASTS1=2501,NASTS2=2501,PI=0.0D0,RHO=2.5D0,MRHO=1.0D0-RHO,MPI=1.0D0-PI)
	PARAMETER (MBAD=0.0D0)
	REAL*8 ASSETS(NASTS2),UTIL1(NASTS1,NASTS1,4,2),V1(NASTS1,4),UTIL2(NASTS2,NASTS2,2,2),V2(NASTS2,2)
	REAL*8 VPR1(NASTS1,4),VPR2(NASTS2,2),SHIRKERS1,UNEMPLOYED1,SHIRKERS2,UNEMPLOYED2
	REAL*8 F1(NASTS1,4),F2(NASTS2,2),FPR1(NASTS1,4),FPR2(NASTS2,2)
	INTEGER DAST(NASTS2,4),DOFFER(NASTS2,4),THETAFLAG
	INTEGER IY,MAXY,YEAR(YY),OLDTHETA2,ITERT
	REAL*8 WAITSVEC(YY),WAITS,MTAU,VTEMP,TOTON,TOTOR,XXX,UIDUR(YY)
	REAL*8 MEANV1,MEANV2,THETA1VEC(YY),ISVEC(YY),UNEMP(YY),UNEMPDUR(YY),WAIT(YY)
	REAL*8 P11,P12,P23,P22,P21,P33,P34,P31,P44,P41,PEE,PEU,PUE,PUU,TOTASS1,TOTASS2
        REAL*8 DIFFMEANV2,OLDMEANV2,WAITERS1,UIS1,ISS1,ATEST(NASTS1),ATEST0(YY),ATEST1(YY)
	REAL*8 BESTTHETA2,BESTMEANV2,BESTTAU2,BESTTOTASS2,BESTSHIRKERS2
	MAXY=YY
	XXX=1
! calibration is read from file calib.dat
	OPEN(2,FILE="calib.dat")
	DO IY=1,YY
	   READ(2,*)YEAR(IY),THETA1VEC(IY),UIDUR(IY),WAIT(IY),ISVEC(IY),UNEMP(IY),UNEMPDUR(IY),WAITSVEC(IY),ATEST0(IY),ATEST1(IY)
	   END DO
! start of computations
! Initialize f 
! we need that right away in the value function iterations.
!	OPEN(2,FILE="v1.dat")
!	OPEN(3,FILE="f1.dat")
245	FORMAT(F18.10)
153	FORMAT(6I7)
!	TOTO=0.0D0
!	DO J=1,NASTS1
!         ASSETS(J)=DBLE(J-1)*INCAST
!	 DO K=1,4
!	      READ(3,245) F1(J,K)
!	      READ(2,245) V1(J,K)
!	 END DO
!	END DO
!	CLOSE(2)
!	CLOSE(3)
	MSIGMA=1.0D0-SIGMA
	MHBAR=(1.0D0-HBAR)**SIGMA
	DO J=1,NASTS2
         ASSETS(J)=DBLE(J-1)*INCAST
	END DO
	DO J=1,NASTS1
! in complex economy, states are assets and job status (1-4)
! 1=offer, 2=unemployed waiting, 3=unemployed with UI, 4=unemployed under income security
         F1(J,1)=1.0D0/DBLE(NASTS1)*(1.0D0-UNEMP(1))
         F1(J,2)=1.0D0/DBLE(NASTS1)*UNEMP(1)/3.0D0
         F1(J,3)=1.0D0/DBLE(NASTS1)*UNEMP(1)/3.0D0
         F1(J,4)=1.0D0/DBLE(NASTS1)*UNEMP(1)/3.0D0
	 DO K=1,4
	    DOFFER(J,K)=1
	    V1(J,K)=-100.0D0
	    VPR1(J,K)=-100.0D0
	 END DO
	END DO
	DO J=1,NASTS2
         F2(J,1)=1.0D0/DBLE(NASTS2)*(1.0D0-UNEMP(1))
         F2(J,2)=1.0D0/DBLE(NASTS2)*UNEMP(1)
	 DO K=1,2
	    V2(J,K)=-100.0D0
	    VPR2(J,K)=-100.0D0
	 END DO
	END DO	
	 TAU1=0.0D0
	 TOTO=0.0D0
	 DO J=1,NASTS1
	   TAU1=TAU1+WAITSVEC(1)*F1(J,2)+THETA1VEC(1)*F1(J,3)+ISVEC(1)*F1(J,4)
	   TOTO=TOTO+(F1(J,1)+WAITSVEC(1)*F1(J,2)+THETA1VEC(1)*F1(J,3)+ISVEC(1)*F1(J,4))
	 END DO
	 TAU1=TAU1/TOTO
	 MTAU=1.0D0-TAU1
! the big iteration starts!
	DO IY=1,MAXY ! YEARS
! calibrate
	   P12=UNEMP(IY)/UNEMPDUR(IY)/(1.0D0-UNEMP(IY))
	   P11=1.0D0-P12
	   P21=1.0D0/UNEMPDUR(IY)
	   P23=(1.0D0-P21)/WAIT(IY)
	   P22=1.0D0-P21-P23
	   P31=1.0D0/UNEMPDUR(IY)
	   P34=(1.0D0-P31)/UIDUR(IY)
	   P33=1.0D0-P31-P34
	   P41=1.0D0/UNEMPDUR(IY)
	   P44=1.0D0-P41
           PEU=UNEMP(IY)/UNEMPDUR(IY)/(1.0D0-UNEMP(IY))
           PEE=1.0D0-PEU
           PUE=1.0D0/UNEMPDUR(IY)
           PUU=1.0D0-PUE
           ATEST1(IY)=ATEST1(IY)*52.0D0
           ATEST0(IY)=ATEST0(IY)*52.0D0
           DO I=1,NASTS1
              IF (ASSETS(I) .LE. ATEST0(IY)) THEN
                 ATEST(I)=1.0D0
              ELSE
                 ATEST(I)=(ATEST1(IY)-ASSETS(I))/(ATEST1(IY)-ATEST0(IY))
              END IF
              IF (ASSETS(I) .GE. ATEST1(IY)) ATEST(I)=0.0D0
           END DO
	   THETA1=THETA1VEC(IY)
	   WAITS=WAITSVEC(IY)
          THETA1=0.303D0
	   IS=ISVEC(IY)
! first compute complex economy
	ITER=0
	ITERF=3
 	ITERV=3
	DO WHILE ((ITERF+ITERV .GT. 4) .AND. (ITER .LT. 100))
	 ITER=ITER+1
! utilities
! util(assets, assets', job status, offer decision)
	 DO J=1,NASTS1
	  DO K=1,NASTS1
           TOTO=MR*ASSETS(J)-ASSETS(K)
	   VTEMP=TOTO+MTAU
	      IF(VTEMP.GT.0.0D0) THEN
	       UTIL1(J,K,1,1)=(VTEMP**MSIGMA)**MRHO/MRHO
	      ELSE
	       UTIL1(J,K,1,1)=BAD*(1.0D0-VTEMP)
	      END IF
	   VTEMP=TOTO+MTAU*WAITS
	      IF(VTEMP.GT.0.0D0) THEN
	       UTIL1(J,K,1,2)=(VTEMP**MSIGMA)**MRHO/MRHO
	       UTIL1(J,K,2,1)=(VTEMP**MSIGMA)**MRHO/MRHO
	      ELSE
	       UTIL1(J,K,1,2)=BAD*(1.0D0-VTEMP)
	       UTIL1(J,K,2,1)=BAD*(1.0D0-VTEMP)
	      END IF
	   VTEMP=TOTO+MTAU*THETA1
	      IF(VTEMP.GT.0.0D0) THEN
	       UTIL1(J,K,3,1)=(VTEMP**MSIGMA)**MRHO/MRHO
	      ELSE
	       UTIL1(J,K,3,1)=BAD*(1.0D0-VTEMP)
	      END IF
	   VTEMP=TOTO+MTAU*IS*ATEST(J)
	      IF(VTEMP.GT.0.0D0) THEN
	       UTIL1(J,K,4,1)=(VTEMP**MSIGMA)**MRHO/MRHO
	      ELSE
	       UTIL1(J,K,4,1)=BAD*(1.0D0-VTEMP)
	      END IF
	  END DO
	 END DO

! value function iterations
! states: (assets, job status) 
22	 PRINT *,'Value function iterations started'
	 VTEST=10.0D0
	 ITERV=1
	 DO WHILE ((CONVERGE/10.0D0 .LE. VTEST) .AND. (ITERV .LT. 200))
	 DO J=1,NASTS1
	    DO K=1,4
	      V1(J,K)=-100000.0D0
	      DAST(J,K)=1
	    END DO
	    DO L=1,NASTS1
! offer, accept it
	     TOTO=UTIL1(J,L,1,1)+BETA*(P11*VPR1(L,1)+P12*VPR1(L,2))
! offer, refuse it
	     TOTOR=UTIL1(J,L,1,2)+BETA*(P21*VPR1(L,1)+P22*VPR1(L,2)+P23*VPR1(L,3))
             IF(TOTO .GT. V1(J,1)) THEN
	       V1(J,1)=TOTO
	       DOFFER(J,1)=1
	       DAST(J,1)=L
	     END IF
             IF(TOTOR .GT. V1(J,1)) THEN
	       V1(J,1)=TOTOR
	       DOFFER(J,1)=2
	       DAST(J,1)=L
	     END IF
	  END DO
! no offer
	  DO L=1,NASTS1
	     TOTO=UTIL1(J,L,2,1)+BETA*(P21*VPR1(L,1)+P22*VPR1(L,2)+P23*VPR1(L,3))
	     IF(TOTO .GT. V1(J,2)) THEN
		V1(J,2)=TOTO
		DAST(J,2)=L
		END IF
	   END DO
	  DO L=1,NASTS1
	     TOTO=UTIL1(J,L,3,1)+BETA*(P31*VPR1(L,1)+P33*VPR1(L,3)+P34*VPR1(L,4))
	     IF(TOTO .GT. V1(J,3)) THEN
		V1(J,3)=TOTO
		DAST(J,3)=L
		END IF
	   END DO
	  DO L=1,NASTS1
	     TOTO=UTIL1(J,L,4,1)+BETA*(P41*VPR1(L,1)+P44*VPR1(L,4))
	     IF(TOTO .GT. V1(J,4)) THEN
		V1(J,4)=TOTO
		DAST(J,4)=L
		END IF
	   END DO
	 END DO
! refresh V
	 VTEST=0.0D0
	 DO J=1,NASTS1
	  DO K=1,4
	   IF (DABS(V1(J,K)-VPR1(J,K)) .GT. VTEST) THEN
	     VTEST=DABS(V1(J,K)-VPR1(J,K))
	   END IF
	   VPR1(J,K)=V1(J,K)
	  END DO
	 END DO
22	 PRINT *,ITER,' IterV',ITERV,' Test',VTEST,' Convergence', CONVERGE,' Tau',TAU1
	 ITERV=ITERV+1
	 END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Distributions
! states: (assets,job status) 
22	 PRINT *,'Distribution function iterations started'
	 FTEST=10.0D0
	 ITERF=1
	 DO WHILE ((FTEST .GT. CONVERGE/100.0D0) .AND. (ITERF .LT. 1000))
	  DO J=1,NASTS1
	   DO K=1,4
	      FPR1(J,K)=0.0D0
	   END DO
	  END DO
	  DO J=1,NASTS1
	      FPR1(DAST(J,1),1)=FPR1(DAST(J,1),1)+P11*F1(J,1)*(2.0D0-DOFFER(J,1))
	      FPR1(DAST(J,1),2)=FPR1(DAST(J,1),2)+P11*F1(J,1)*(DOFFER(J,1)-1.0D0)
	      FPR1(DAST(J,1),2)=FPR1(DAST(J,1),2)+P12*F1(J,1)
	      FPR1(DAST(J,2),1)=FPR1(DAST(J,2),1)+P21*F1(J,2)
	      FPR1(DAST(J,2),2)=FPR1(DAST(J,2),2)+P22*F1(J,2)
	      FPR1(DAST(J,2),3)=FPR1(DAST(J,2),3)+P23*F1(J,2)
	      FPR1(DAST(J,3),1)=FPR1(DAST(J,3),1)+P31*F1(J,3)
	      FPR1(DAST(J,3),3)=FPR1(DAST(J,3),3)+P33*F1(J,3)
	      FPR1(DAST(J,3),4)=FPR1(DAST(J,3),4)+P34*F1(J,3)
	      FPR1(DAST(J,4),1)=FPR1(DAST(J,4),1)+P41*F1(J,4)
	      FPR1(DAST(J,4),4)=FPR1(DAST(J,4),4)+P44*F1(J,4)
	  END DO
	  FTEST=0.0D0
	  TOTO=0.0D0
	  DO J=1,NASTS1
	   DO K=1,4
	     IF (DABS(F1(J,K)-FPR1(J,K)) .GT. FTEST) THEN
	         FTEST=DABS(F1(J,K)-FPR1(J,K))
	        END IF
	        F1(J,K)=FPR1(J,K)
	        TOTO=TOTO+F1(J,K)
	   END DO
	  END DO
22	  PRINT *,ITER,' IterF',ITERF,' Test',FTEST,' Convergence',CONVERGE/100.0D0,
22     1     ' Total F',TOTO
	  ITERF=ITERF+1
	 END DO !F iteration
          TUTU=0.0D0
          DO J=NASTS1-5,NASTS1
             DO K=1,4
                TUTU=TUTU+F1(J,K)
             END DO
          END DO
          IF(TUTU .GT. 0.01D0) THEN
             PRINT *,'Grid border attained by F1',TUTU
          END IF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! new tax rate, given shirkers
	 TOTO=0.0D0
	 TOTON=TAU1
	 TAU1=0.0D0
	  SHIRKERS1=0.0D0
	  UNEMPLOYED1=0.0D0
	  WAITERS1=0.0D0
	  UIS1=0.0D0
	  ISS1=0.0D0
	   DO J=1,NASTS1
	     DO K=1,4
	      IF(DOFFER(J,K) .EQ. 2) THEN
	      SHIRKERS1=SHIRKERS1+F1(J,K)
	      END IF 
	     END DO
	     WAITERS1=WAITERS1+F1(J,2)
	     UIS1=UIS1+F1(J,3)
	     ISS1=ISS1+F1(J,4)
	     UNEMPLOYED1=UNEMPLOYED1+F1(J,4)+F1(J,2)+F1(J,3)
	     TAU1=TAU1+THETA1*F1(J,3)+IS*F1(J,4)*ATEST(J)+WAITS*F1(J,2)
	     TOTO=TOTO+F1(J,1)
	    END DO
	 TAU1=TAU1/(TOTO+TAU1)
	 MTAU=1.0D0-TAU1
	END DO !BIG iteration
	TOTASS1=0.0D0
	MEANV1=0.0D0
	DO J=1,NASTS1
	  DO K=1,4
	    MEANV1=MEANV1+F1(J,K)*V1(J,K)
	    TOTASS1=TOTASS1+F1(J,K)*ASSETS(J)
	    IF ((F1(J,K) .GT. 0.001D0) .AND. (V1(J,K) .LT. -1000.0D0)) THEN
	       PRINT *,J,' ',K,' ',F1(J,K),' ',V1(J,K)
	       END IF
	 END DO
	END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! and now compute simple economy
! initialize using complex economy for first pass
	IF (IY .EQ. 1) THEN
	OPEN(4,FILE="results2.test")
100	FORMAT('year pi theta1 theta2 tau1 tau2 meanv1 meanv2 totass1
     C   totass2 shirker1 shirkers2')
	WRITE(4,101)
101	FORMAT('============================================')
	DO J=1,NASTS1
	   F2(J,1)=F1(J,1)
	   F2(J,2)=F1(J,2)+F1(J,3)+F1(J,4)
	   VPR2(J,1)=V1(J,1)
	   VPR2(J,2)=V1(J,3)
	   END DO
	DO J=NASTS1+1,NASTS2
	   F2(J,1)=0.0D0
	   F2(J,2)=0.0D0
	   VPR2(J,1)=V1(NASTS1,1)
	   VPR2(J,2)=V1(NASTS1,3)
	   END DO
	END IF
	THETA2=THETA1
        OLDTHETA2=THETA2
	TAU2=0.0D0
	TOTO=0.0D0
	DO J=1,NASTS2
	     TAU2=TAU2+THETA2*F2(J,2)
	     TOTO=TOTO+F2(J,1)
	    END DO
	 TAU2=TAU2/(TOTO+TAU2)
        DIFFMEANV2=10.0D0*DABS(MEANV1)
	MTAU=1.0D0-TAU2
	MEANV2=0.0D0
	VTEST=0.0D0
	ITERT=0
	THETAFLAG=0
        THETA2=THETA2+0.01D0
        DO WHILE ((THETAFLAG .EQ. 0) .AND. (ITERT .LT. 1000) .AND. (THETA2 .LT. 1.0D0) .AND. (THETA2 .GT. -0.201D0))
           THETA2=THETA2-0.01D0
           OLDMEANV2=MEANV2
	   ITERT=ITERT+1
	ITER=0
	ITERF=3
 	ITERV=3
	DO WHILE ((ITERF+ITERV .GT. 4) .AND. (ITER .LT. 100))
	 ITER=ITER+1
! utilities
! util(assets, assets', job status, offer decision)
	 DO J=1,NASTS2
	  DO K=1,NASTS2
           TOTO=MR*ASSETS(J)-ASSETS(K)
	   VTEMP=TOTO+MTAU
	      IF(VTEMP.GT.0.0D0) THEN
	       UTIL1(J,K,1,1)=(VTEMP**MSIGMA)**MRHO/MRHO
	      ELSE
	       UTIL1(J,K,1,1)=BAD
	      END IF
              VTEMP=TOTO+MTAU*IS*ATEST(J)
              IF(VTEMP.GT.0.0D0) THEN
               UTIL1(J,K,1,2)=(VTEMP**MSIGMA)**MRHO/MRHO
              ELSE
               UTIL1(J,K,1,2)=BAD
              END IF
	   VTEMP=TOTO+MTAU*THETA2
	      IF(VTEMP.GT.0.0D0) THEN
	       UTIL1(J,K,3,1)=(VTEMP**MSIGMA)**MRHO/MRHO
	      ELSE
	       UTIL1(J,K,3,1)=BAD
	      END IF
	  END DO
	 END DO

! value function iterations
! states: (assets, job status) 
22	 PRINT *,'Value function iterations started'
	 VTEST=10.0D0
	 ITERV=1
	 DO WHILE ((CONVERGE .LE. VTEST) .AND. (ITERV .LT. 200))
	 DO J=1,NASTS2
	    V2(J,1)=-100000.0D0
	    V2(J,2)=-100000.0D0
	    DO L=1,NASTS2
! offer, accept it
	     TOTO=UTIL1(J,L,1,1)+BETA*(PEE*VPR2(L,1)+PEU*VPR2(L,2))
! offer, refuse it
	     TOTOR=UTIL1(J,L,1,2)+BETA*(PUE*VPR2(L,1)+PUU*VPR2(L,2))+BAD
             IF(TOTO .GT. V2(J,1)) THEN
	       V2(J,1)=TOTO
	       DOFFER(J,1)=1
	       DAST(J,1)=L
	     END IF
             IF(TOTOR .GT. V2(J,1)) THEN
	       V2(J,1)=TOTOR
	       DOFFER(J,1)=2
	       DAST(J,1)=L
	     END IF
	  END DO
! no offer
	  DO L=1,NASTS2
	     TOTO=UTIL1(J,L,3,1)+BETA*(PUE*VPR2(L,1)+PUU*VPR2(L,2))
	     IF(TOTO .GT. V2(J,2)) THEN
		V2(J,2)=TOTO
		DAST(J,2)=L
		END IF
	   END DO
	 END DO
! refresh V
	 VTEST=0.0D0
	 DO J=1,NASTS2
	  DO K=1,2
	   IF (DABS(V2(J,K)-VPR2(J,K)) .GT. VTEST) THEN
	     VTEST=DABS(V2(J,K)-VPR2(J,K))
	   END IF
	   VPR2(J,K)=V2(J,K)
	  END DO
	 END DO
22	 PRINT *,ITER,' IterV',ITERV,' Test',VTEST,' Convergence', CONVERGE,' Tau',TAU2
	 ITERV=ITERV+1
	 END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Distributions
! states: (assets,job status) 
22	 PRINT *,'Distribution function iterations started'
	 FTEST=10.0D0
	 ITERF=1
	 DO WHILE ((FTEST .GT. CONVERGE/100.0D0) .AND. (ITERF .LT. 1000))
	  DO J=1,NASTS2
	   DO K=1,2
	      FPR2(J,K)=0.0D0
	   END DO
	  END DO
	  DO J=1,NASTS2
	      FPR2(DAST(J,1),1)=FPR2(DAST(J,1),1)+PEE*F2(J,1)
	      FPR2(DAST(J,2),2)=FPR2(DAST(J,2),2)+PUU*F2(J,2)
	      FPR2(DAST(J,2),1)=FPR2(DAST(J,2),1)+PUE*F2(J,2)
	      FPR2(DAST(J,1),2)=FPR2(DAST(J,1),2)+PEU*F2(J,1)
	  END DO
	  FTEST=0.0D0
	  TOTO=0.0D0
	  DO J=1,NASTS2
	   DO K=1,2
	     IF (DABS(F2(J,K)-FPR2(J,K)) .GT. FTEST) THEN
	         FTEST=DABS(F2(J,K)-FPR2(J,K))
	        END IF
	        F2(J,K)=FPR2(J,K)
	        TOTO=TOTO+F2(J,K)
	   END DO
	  END DO
22	  PRINT *,ITER,' IterF',ITERF,' Test',FTEST,' Convergence',CONVERGE/100.0D0,
22     1     ' Total F',TOTO
	  ITERF=ITERF+1
          TUTU=0.0D0
          DO J=NASTS2-5,NASTS2
             DO K=1,2
                TUTU=TUTU+F2(J,K)
             END DO
          END DO
          IF(TUTU .GT. 0.01D0) THEN
             PRINT *,'Grid border attained by F2',TUTU
          END IF
	 END DO !F iteration
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! new tax rate, given shirkers
	 TOTO=0.0D0
	 TOTON=TAU2
	 TAU2=0.0D0
	  SHIRKERS2=0.0D0
	  UNEMPLOYED2=0.0D0
	   DO J=1,NASTS2
	     DO K=1,2
	      IF(DOFFER(J,K) .EQ. 2) THEN
	      SHIRKERS2=SHIRKERS2+F2(J,K)
	      END IF
	     END DO
	     UNEMPLOYED2=UNEMPLOYED2+F2(J,2)
	     TAU2=TAU2+THETA2*F2(J,2)
	     TOTO=TOTO+F2(J,1)
	    END DO
	 TAU2=TAU2/(TOTO+TAU2)
!	 TAU2=(TAU2+TOTON)/2.0D0
	 MTAU=1.0D0-TAU2
	END DO !BIG iteration
	TOTASS2=0.0D0
	MEANV2=0.0D0
	DO J=1,NASTS2
	  DO K=1,2
	    MEANV2=MEANV2+F2(J,K)*V2(J,K)
	    TOTASS2=TOTASS2+F2(J,K)*ASSETS(J)
	 END DO
	END DO
	PRINT *,ITERT,' Year:',YEAR(IY),' Mean V1:',MEANV1,' Mean V2:',MEANV2,' Theta1:',THETA1,' Theta2:',THETA2, ' Tau1:',TAU1,' Tau2:',TAU2,' Shirkers1:',SHIRKERS1,' Unemployed1:',UNEMPLOYED1,"=",WAITERS1,UIS1,ISS1
	      OLDTHETA2=THETA2
        IF (DABS(MEANV1-MEANV2) .LT. DIFFMEANV2) THEN
           DIFFMEANV2=DABS(MEANV1-MEANV2)
	   BESTTHETA2=THETA2
           BESTMEANV2=MEANV2
           BESTTAU2=TAU2
           BESTTOTASS2=TOTASS2
           BESTSHIRKERS2=SHIRKERS2
        ELSE
           IF (DABS(MEANV1-MEANV2) .GT. 5.0D0*DIFFMEANV2) THEN
           THETAFLAG=1
           END IF
        END IF
        PRINT *,ITERT,' Year:',YEAR(IY),' Mean V1:',MEANV1,' Mean V2:',MEANV2,' Theta1:',THETA1,' Theta2:',THETA2,' IterV:',ITERV
	END DO !END ITERATION OVER THETA2
	WRITE(4,102) YEAR(IY),PI,THETA1,BESTTHETA2,TAU1,BESTTAU2,MEANV1,BESTMEANV2,TOTASS1,BESTTOTASS2,SHIRKERS1,BESTSHIRKERS2
 102	FORMAT(I4,3F6.3,2F9.6,2F12.5,2F9.5,2F6.3)
	END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! save v and f matrices
	OPEN(2,FILE="v1.test")
	OPEN(3,FILE="f1.test")
	DO J=1,NASTS1
	 DO K=1,4
	      WRITE(3,245) F1(J,K)
	      WRITE(2,245) V1(J,K)
          END DO
	END DO
	OPEN(5,FILE="v2.test")
	OPEN(6,FILE="f2.test")
	DO J=1,NASTS2
	 DO K=1,2
	      WRITE(5,245) F2(J,K)
	      WRITE(6,245) V2(J,K)
          END DO
	END DO
	PRINT *,'All done!' 
	END
